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ABSTRACT 



The object of this thesis is to examine, via simulation, 
the properties of a class of single-server, first-come 
first-served queues in which the service times and inter- 
arrival times, while still marginally exponential, are 
auto-correlated and cross-correlated. The correlation is 
introduced by letting the service and interarrival sequences 
be EARMA-type processes, where EARMA stands for exponential 
autoregressive, moving average sequence. An extension of 
these ideas brings in the cross-correlation. The waiting 
times in the dependent queue are compared to the waiting 
times (with known distribution) in the independent M/M/1 
queue using data analytic and formal statistical methods. 
Variance reduction techniques are also studied; these use 
distributionally known aspects of the M/M/1 queue (simulated 
with the same exponential sequences as the dependent queue) 
to control the unknown aspects of the dependent queue. 
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I. 



INTRODUCTION 



The object of this thesis is to examine, via simulation, 
the properties of a class of single-server first-come first- 
served queues in which the service times and interarrival 
times, while still marginally exponential, are auto-correlated 
and cross-correlated. The correlation is obtained by making 
the service and interarrival time sequences multivariate EARMA 
processes. An ancillary aspect of the study is the investi- 
gation of the efficiency of a number of control variable 
variance reduction schemes. These schemes are made possible 
by the simplicity of the structure of the multivariate EARMA 
processes, and the fact that the null case (no correlation) 
gives the M/M/1 cjueue , whose properties are known analytically. 

An immediate problem which arises in the simulation is 
that it is not known how far out along the sample path the 
simulation must be carried in order for the steady state 
properties of the queue to hold. This is a general problem 
which arises in queuing simulations and, in order to cope with 
it, the simulations are run in 500 parallel replications. 

This allows one to use many methods of time series analysis 
and data analysis to obtain graphical verification in a se- 
quential manner of the approach to a steady state. The repli- 
cations also allow for examination of the complete steady 
state waiting time distribution in the queue , and not merely 
the mean waiting time. 
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While all of these aspects of the queuing simulation 
are interwoven, we can separate them out for purposes of 
discussion as follows. 

(i) An elementary aspect of the queue which can be examined 
is the mean waiting time. The most efficient (smallest 
variance) estimate of this is the sample path average 
which is again averaged over replications. This repli- 
cation of course reduces the variance by a factor of 

m, where m is the number of replications. Furthermore 
the existence of replications allows one to look at the 
complete distribution of the sample path average esti- 
mate. This distribution should be normally distributed, 
possibly even before a steady state has been reached. 

The object of this part of the study is to see how 
the mean steady state waiting time in the queue varies 
with the traffic intensity, t, and the correlation 
parameters, and in particular, how this mean waiting 
time differs from the waiting time in the M/M/1 queue 
(null case) . 

(ii) A more detailed aspect of the queue which has to be 
examined is the complete distribution of the steady 
state waiting time. In the M/M/1 queue this is 
exponential, given that the waiting time is greater 
than zero. 

The first problem here is to determine that one 
is actually looking at the steady state waiting time. 
This is done by examining the empirical distribution 
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function of the waiting times across replications 

at times n 1 and n 2 (n 2 > n^) , where the index n 

refers to the customer number. Of course the two 

samples from W and W , each of size m, will be 
n l n 2 

correlated and this makes standard statistical methods 
for comparing samples invalid. However the correla- 
tion between these samples can be calculated in order 
to determine that n 2 is sufficiently large relative to 
n^ for the correlation to be negligible. In actual 
fact the comparison is done by looking at box plots 

of the W ' s at successive values n. < n_ < n_ < ... . 
n 12 3 

This method is graphic and allows the simulator to 
see how the steady state is being approached. 

(iii) A third problem which is involved in looking at the 
complete distribution of the waiting time is that a 
sample of size m = 500 is not sufficient to really 
get a good estimate of the steady state distribution 
of the waiting time. Thus successive samples along 
the parallel sample paths are combined, where it is 
determined that the samples are far enough apart to 
be approximately independent. 

Again the object here is to see how the complete 
distribution of the steady state waiting time differs 
from the exponential distribution (null M/M/ 1 case) 
as traffic intensity, t, and correlation parameters 
vary. In particular heavy traffic theory says that 
this distribution should be approximately exponential 
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as the traffic intensity parameter approaches 1. 

The difficulty in verifying this, and the heavy 
traffic approximation to the mean waiting time, is 
that the time to reach the steady state becomes very 
long as t + 1. An additional question that arises 
then is whether it is possible to start the simulation 
in such a way (approximate steady state conditions) 
so as to speed up this convergence. It should be 
remembered, that the queue with correlation is not 
regenerative and no initial conditions for stationarity 
are known or are likely to be found. 

(iv) Variance reduction is always helpful in a simulation 
of this kind since real detail is always masked by 
sampling fluctuations i.e., the variability in the 
estimate of the unknown mean waiting time. Because 
of the simple probabilistic-linear structure of the 
EARMA processes, running on M/M/1 queue in parallel 
with the correlated queue using common exponential 
variates provides a very good control variable for the 
waiting time if the correlation, in a rough sense, 
in the EARMA queue, is small. However, as this 
correlation increases, the amount of control decreases. 
Thus several other control variables are examined and 
combined into a multiple control for the correlated 
queue. In this way it is hoped to obtain variance 
reduction for any values of the correlation parameters 
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in the EARMA queue. The actual degree of control 
which can be attained is examined, using the repli- 
cated simulation, as a function of the parameters 
of the queue. 
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II. THE EXPONENTIAL AUTOREGRESSIVE-MOVING 
AVERAGE PROCESS , EARMA (1,1) 



The first order moving average process and the first- 
order autoregressive processes were combined by Jacobs and 
Lewis (1977) to form the EARMA (1,1) sequence. We first 
describe these two first order sequences, and then the method 
of combining them. 

A. THE EXPONENTIAL AUTOREGRESSIVE PROCESS, EAR ( 1) 

The standard linear, first-order autoregressive model for 
a stationary sequence of random variables {X^} is defined 
by the equation 



where p is a constant which is less than 1 in absolute value 
and {e^} is a sequence of independent and identically dis- 
tributed random variables. Gaver and Lewis (1978) showed 
that if the {X^} sequence were to have an exponential marginal 
distribution with parameter a , then the parameter p should be 
greater than or equal to zero and less than one, and 
should be zero with probability p or an exponential (A) 
random variables, E^, with probability 1-p. Thus 



(II. 1) 




pX i_l + ; i = 0, ±1, ±2 




pX i _ 1 + e i ; i = 0, ±1, ±2 



becomes 
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with probability p 



(11. 2) X i = 

pX^_^ + E^ with probability 1-p, 

i = 0, ±1, ±2, ... where {E^} is an identical independent 

distribution (i.i.d.) sequence of exponential (A) random 
variables. Note that for this EAR(l) model the distribution 
of the depend on p, the multiplicative weight of X^_^. 

Also is not an absolutely continuous random variable. 

Thus standard results (Mallows, 1968) to the effect that the 
X^ in an autoregressive process become approximately normal 
as p 1 do not hold. In fact in the EAR(l) process (II.2), 
the X^ are exponentially distributed for 0 <_ p < 1. 

B. THE EXPONENTIAL MOVING AVERAGE PROCESS, EMA(l) 

The first order moving average exponential process EMA(l) 
was developed by Lawrance and Lewis (1977) in the form of 

6E i w.p. 8 

(11. 3) X. = 

l 

8E i + E^_ x W.p. (1-6) , 

i = 0, ±1, ±2, ... , where 6 is greater than or equal to zero 

and less than or equal to one, and {E^} is again a sequence 
of i.i.d. exponential (A) random variables. (This is the 
backward case — a forward case is also possible.) The X^'s 
have an exponential marginal distribution and are only serially 
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dependent for lag one; this model is highly tractable and 
a full account of the statistically useful properties was 
obtained by Lawrance and Lewis (1977) . The forward model 
combines E^ with E^ + ^ instead of with E^_^. 

C. THE EXPONENTIAL AUTOREGRESSIVE-MOVING AVERAGE 
PROCESSES EARMA (1,1) 

The EAR (1) model of Section IIA and the EMA(l) model 
of Section IIB are combined in the following form to give 
an EARMA (1,1) process. This is a non-Markovian model which 
contains the EMA(l) and EAR (1) models as a special case. The 
defining equations for the EARMA (1,1) process are 

$E^ , w.p. 3, 

(II. 4) X j _ = 

3Ef + A i _ 1 , w.p. 1-3, 

(0 < 3 < 1; i = 0, ±1, ±2 , . . . ) 



where 



pA i-2 w.p. p, 

A i-1 = 

pA i-2 + E i _ 1 / w.p. 1-p, 

(0 < p < 1; i = 0, ±1, ±2, ...) 

and {E^ } is, as usual, a sequence of i.i.d. exponential (X) 
random variables. Thus instead of E^ being combined with the 
exponential random variable E^_^, it is combined with an 
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exponentially distributed random variable which is an EAR (1) 
combination of E^_^, E j__ 2 ' E i- 3 ' ••• • The serial correla- 

tions for the EARMA (1,1) process are given in Jacobs and 
Lewis (1977) as 

p(j) = corr (X^X^ + j ) = c(p/B) p ^ 1 

(j = 0,11,12,13, 0 _< p < 1; 0 _< 3 _< 1). 

where c(p,3) = 3(1-6) (1“P) + (l-3) 2 p. 

Note that when p = 0 the EARMA (1,1) process reduces to an 
EMA(l) process; if 3 = 0 it is the autoregressive EAR(l) 
process, and when 3 = 1 it is the usual sequence of 
exponential (A) random variables. We will say that 
is autoregressive over E^ , E^_^, ... . 
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III. USE OF MIXED STRUCTURES EARMA (1,1) 
IN MODELLING QUEUES 



Consider for simplicity a queue with a single input 
stream and a single server and a f irst-come-f irst-served 
(FIFO) service discipline. Let S^, i=0, 1, 2, ... , 
denote the service time for the i th arrival, and let X^, 
i = 1, 2, ... , denote times between arrival of the ith and 

(i-l)th customers (interarrival times). As is usual we assume 
that the first customer (with service time S Q ) arrives at 
time zero and finds the queue empty. 

If the {S^} and {X^} sequences are i.i.d. exponential 
random variables with parameters A and a respectively, we 
have M/M/1 queue. 

Now let 

E^ be i.i.d. exponential (A); i= 0, 1, 2, ... , 
be i.i.d. exponential (a); i=0, 1, 2, ... . 

We want to model queues with correlated (autocorrelated 
and/or cross-correlated) service and inter-arrival times, 
the service and inter-arrival times both having marginally 
exponential distribution (MD^4D/1 queue). Exponential 
marginals is a common assumption and with it the MDxMD/1 
queuing model includes the M/M/1 queue as a special case. 
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The dependence in the arrival and service processes 
is created here by defining and as a bivariate depen- 
dent sequence of random variables with exponential marginal 
distributions. This is done by letting {S^} be EARMA (1,1) 
over , (a/A)e^, (a/A)e^ , ... , i.e. 



S 



0 



E 



0 



6E X w.p. B; A^ = (a/A)e 1 ; 

6E l + A 1 , w.p. (1-6) ; 



BE, 



w.p. 3; 



pA 1 



w.p. p ; 



‘ A 2 = 



BE 2 + A 2 , w.p. (1-6) ; 



a 



pA x + jz 2 w.p. (1-p) ; 



BE i , w.p. 6; 

s . = 

1 



pA i-l W ‘ P * p; 

A. = 

1 



BEi + A i , w.p. (1-6) 7 



pA 



i-1 




w.p . 



(1-p) 



We also let 



X. = 
i 



i = 1, 2, 3, 



/ 
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although one could further make the X^ ' s correlated. With 
the present assumptions the input process is still a Poisson 
process, as in the M/M/1 queue. The cross-coupling between 
the sequence {S^} and {X^} is apparent from the structure. 

Note that the S^'s are now correlated because of the cross- 
coupling. The {S^} is a sequence of dependent exponential 
random variables, but not an EARMA(1,1) sequence. It is a 
type of pseudo-EARMA (1, 1) sequence. 

The interpretation of this queuing model is that the 

server tends to speed up if the queue gets long (in the 

\ 

past) . Of course he also slows down when the queue gets 
short. In the case where p = 0 and 8=0 then the correlation 
between {S^} and {X^} will be equal to one, i.e., p 
S i = (a/X)X^. Also when 8=1/ and for any value of p, 
the queuing model will be identical to the M/M/1 queue. Thus 
the M/M/1 is included as a special case. Some more complicated 
schemes for coupling the service and inter-arrival processes 
are given in Lewis and Shedler (1978) . Jacobs (1978a, 1978b) 
has studied heavy traffic approximations for EARMA-type queues. 

Other attempts to model M/M/1 (FIFO) type queues with 
correlated service and interarrival times have been made. 

None of them, however, have given fixed (exponential) 
marginal distributions for and X^, or cross-coupling between 
these sequences. There are, of course, Markovian schemes 
based on birth-death representations which give FIFO queues 
in which the service times and interarrival times are correlated 
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(Cox and Smith, 1961) , but the scheme is quite different 
from that given here. In particular the interarrival and 
service times are not exponentially distributed, and this 
may not conform to what is observed in practice. 
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IV. THE SIMULATION MODEL 



A. INTRODUCTION AND NOTATION 

Since both queues we are considering, the M/M/1 queue 
and the MDxMD/1, are single-server, first-in-first out queues 
the waiting times W n in both are given by the recursion 
equation 



W = 0 

o 

W , , = max (W + S - X ,, ,0) ; n = 0,1,2, .. 

n+1 n n n+1 

These equations are used to generate successive W n 1 s in the 
simulation. This is the only aspect of the queue which will 
be considered. To have looked at quantities like the number 
of customers in the queue at departure times would have com- 
plicated the programming and taken away from the primary aim 
of the thesis, which is to explore the effect of correlation 
on the queue and to use some relatively modern statistical 
methodology to do this. 

We denote the waiting time of the nth customer to 
arrive in the jth realization of the queue by W (j). The 
usual sample path estimate for the mean, E(w), of the 
limiting distribution of waiting times, F (x) , where 

F (x) = lim P(w < x) , 

w n — 
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is tile sample path average 



n 



W (j) 
n J 



i i 

n L 

1=0 



u > 



which can be shown to converge to the mean E(W) as n -*■ ®>. 

The M independent sample path realizations W (j), 
n= 0,1,2, ..., j = 1,2, ..., M can be used to obtain estimates 

of 

(1) the distribution of the estimate W (j). This should 

n J 

be normal for large n. This can be examined from 
the sample W (j), j = 1, In addition the 

mean of this sample 

M M n 

W =ilw(j) = j 7 W. (j) 

n M L n J Mn L L l J 

j = l j = l Jl=0 

is a grand mean estimating E(w) whose variance is 
computable from the sample W^(j), j = 1, ..., M. 

Thus 

Var(W^) = ^(sample variance of W^(j)'s) 

M 

= H<MTT ^ * V 2) 

j = l 
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( 2 ) 



The distribution of W (not only its mean) from the 

n 

sample W n (j), j = 1, . . . , M. Since this is a random 

sample, all the classical methods are available, 

e.g., histograms, normal plots, empirical c.d.f's. 

(3) The correlations between successive waiting times, 

such as corr (W (j),W ( j ) ) , n. < n„ . This is 

n^ n.2 12 

given as 



1 

M 



M 

( I (W 

j = l 



n. 



(j) 



) (W„ (j ) - V7 ) ) 



n 



n 



n. 



divided by the sample standard deviations of the two 
samples, where 



W 



n 



1 



1 

M 



M 

I (W (j)) 
, n l 



is an average across replications. 

B. PROGRAM STRUCTURE 
1 . General 

To simulate this particular EARMA (1,1) with cross- 
correlated service time and inter-arrival time, a FORTRAN 
program has been written to calculate the mean waiting time 
W n (j). It also calculates the mean of various statistics 
for the M/M/1 queue using the same exponential deviates that 
are used to generate the correlated queue. These statistics 
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are used to control the estimates from the MDxMD/1 queue for 
purposes of variance reduction. The maximum number of 
customers which the program can handle for each replication 
is n = 10,000. This computation is divided into a maximum of 
10 steps and the samples of waiting times at these points are 
stored in a file. Thus for example, if we use the maximum 
size n = 10,000, we have ten data sets W ]_qoo ^ ' W 2000 ^' 

..., W 1Q 000 (j), j = 1, ..., M. Thus the time evolution of 
the queue can be studied. Moreover there is a facility to 
restart the simulation to get Q00 (j), W 12 000 ^ ••• 

This way the evolution of the queuing process can be studied 
as far out as needed, e.g., to n = 270,000, or n = 1,000,000 etc. 
Another program will take care of reading from the output file 
to do data analysis, i.e., box .plot of waiting times at 
various steps, average waiting times W N (j) and their statistics. 

2 . MAINl Program 

a) The variables are to be read in as follows: 



N 



RS, RS 
BETA, RHOS 
KN 
KS 



NREP 

NSTEP 



= Number of arrivals to be generated 
(max is 10,000) 

= Arrival, Service rate. 

= 3, P 

- Number of the run 

= 1 if first run, so that n = 0,1,..., N 

0 otherwise (i.e., restart) 

= Number of replications (max 500) 

= Number of steps in this run 
(maximum of 10) 
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NPAST = Number of arrivals of last run 

(i.e., last run stopped at n = 120,000) 

IEX, IES = Seeds to generate exponential random 

variable of Arrival, Service times 

IBRTA, IRHOS = Seeds to generate uniform random 

variables to generate probabilities 
in the EARM(1,1) process of moving 
average and auto-regressive. 

b) The output will be stored in two separated files. 

i) For continuation of the generation of the waiting 

times of each replication to the next run (i.e. , the restart) , 
the variables stored are as follow: 



Correlated queue 



ARIV 




cumulative inter-arrival time (these are 
the same for both queues) 


SERV 


= 


cumulative service time 


SNEW 


= 


service time of the last customer in this 
run 


WAIT 


= 


waiting time of the last customer in this 
run 


WBAR 


= 


cumulative waiting time 


ALAST 


= 


Auto-regressive component of the inter- 
arrival time of the last customer in this 
run 


TL 




Number of customers who have arrived since 
the last zero waiting time, i.e., backward 
recurrence time in the renewal process of 
times of emptiness in the queue 


TR 


= 


Number of zero waiting times up to the 
termination point of this run. 



Uncorrelated queue 



The variables are defined as the same as for 
the correlated queue; in the program the quantities are 
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suffixed by M, e.g. WAITM is the waiting time for the 
uncorrelated queue and WAIT the waiting time for the 
correlated queue. 

i-i) For analyses and displays of the data in 
each run and each step the data, including the data with 
control variables added, are stored in another file in the 
sequence of waiting time in each replication of the corre- 
lated queue, average waiting time in each replication of 
the correlated queue, waiting time in each replication of the 
uncorrelated queue, average waiting time in each replication 
of M/M/1 queue, waiting time (controlled) and average waiting 
time (controlled) of the correlated queue. 

' 3 . Subroutine CARMA 

This subroutine is used to create a string of service 
times, S n , which are cross-correlated to the inter-arrival 
time X n , and also to create the moving average structure. 

4. Subroutine CONVAR, COPYSM, COPYMV and MPROCV 
These four subroutines are used to compute the 

statistics from the M/M/1 queue which are used to control 

the data from MDxMD/1 queue for purposes of variance reduction. 

5 . MAIN2 Program with Subroutine COMPAR, FILL and STAT 
These programs and subroutines are used to display 

and analyze the data for each run. They produce the figures 
given in this thesis. (See Appendix II.) 
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a. Description of Box Plot 

Subroutine COMPAR will produce vertical box 
plots (McNeil 1977) parallel to each other and write out 
the minimum and the maximum value of all the data in the 
array. 

Each box plot is obtained by first calculating 
the lower and upper quartiles and the median of the batch 
of numbers. 

Where the median of a batch of numbers is the 
value for which half of the numbers in the batch are larger 
and half are smaller, the lower quartile is the value that 
divides the batch into two parts, with 1/4 of the numbers 
below this value, and 3/4 above it. Similarly the upper 
quartile is the value for which 3/4 of the numbers are 
below it and 1/4 above it. 

The narrow rectangular box with ends corres- 
ponding to the lower and upper quartiles contains the median 
point (*) . The height of the box, the interquartile 
distance , D , is then measured out on each side of the 
quartiles. Then the lowest and highest data points that 
fall between these quantities are also marked by crosses (X) . 
Finally, any numbers whose positions are outside these crosses 
are marked with circles (0), those more than 1.5 inter- 
quartile distances outside each quartile getting the 
symbol (@) . The digit besides the "outlier" positions is, 
to indicate the number of data points that are too close 
to be marked separately. (See Figure 1 . ) 
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FIGURE 1. Description of Box Plot 



V. ANALYSIS OF SIMULATION RESULTS 



A. INTRODUCTION 

The object of the simulation was to examine the effect 
of the cross correlation, indexed by the coefficients p and 
8 in the EARMA (1,1) process, on the stationary waiting time 
W for various values of the traffic intensity t. Both the 
mean and distribution of W were compared to the known mean 
and distribution of the waiting time in an M/M/1 queue 
(same t and same mean service time) . 

To perform this comparison we initially simulated this 
model for values .25, .50, .75, .90 of the taffic intensity 
(t) , for values 0.0, .25, .50, .75, & .90 of the autoregressive 

parameter (p) , and the single value 0.50 for the moving aver- 
age coefficient (8) for each combination of t and p. These 
values were chosen to limit the amount of computing to a 
physically feasible range, with the idea that higher values 
of p and t would be examined if time permitted. The number 
of customers, n, until steady state is reached goes up drama- 
tically as p and t increase. The simulations are run in 500 
parallel replications, and the 10,000 arrivals for each run 
are divided into 10 steps of 1,000 arrivals. 

Successive groups of 10,000 arrivals were run for each 
combination of p and t until it was determined from an analy- 
sis of the output that the simulation had reached a steady 
state. This determination was made from box-plots of samples 
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M, for n 2 > n^, from plots of 



W_ ( j ) , W ( j ) , j = 1 , 

.1 _ n 2 

W n and W n# . and from formal statistical comparisons of the 
two samples with n 2 large enough compared to n^ that the 
correlation between the two samples is negligible. 

B. ESTIMATION OF TIME TO STEADY STATE 

For each run, the distribution of the sample path average 
estimate W (j) for each step is observed by the method of 
box plots. By comparing these plots along the n axis at 
steps of arrivals of 1,000, convergence to steady-state and 
normality can be determined. Figures la-g show the conver- 
gence to steady state of the M/M/1 queue at t = .75. For 
the EARMA-type MDxMD/1 queue with t= .75, p = .75, 6 = .50, 
the box plots are shown in Fig. 2a-g. The number of arrivals 
at the estimated steady state for the M/M/1 queue and for 
the correlated queue for any combination of p and t are 
shown in table 1. 

In Figures 1 and 2 note that the variance of W (j) is 
decreasing as 1/n, so that there is a continual shrinkage 
from left to right in the plots. The symmetry and lack of 
extreme values gives an indication of normality of the esti- 
mates. Note that the scales on successive figures change 
through la to lg etc. This is because the plots are not 
commensurate. Range of values is indicated by the values 
at top and bottom on the left, e.g., 2.3460522 and 0.3639930 
in Figure la. 
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Traffic intensity, t 





.25 


. 50 


.75 


.90 


.00 


40,000 


40,000 


50,000 


50,000 


.25 


40,000 


50,000 


50,000 


70,000 


. 50 


40,000 


50,000 


50,000 


70,000 


.75 


40,000 


50,000 


70,000 


100,000 


.90 


60,000 


70,000 


100,000 


100,000 



Table 1. Number of arrivals that the simulation 
run indicates are needed approximately 
to reach steady state (for 8=0.5). These 
are read off of plots such as those given 
for t = 0.75, p = 0.75 and 8 = 0.5 in 
Figures 2a-g. 
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Figure la. Box plots derived from 500 replications, j=l,...,500 of the 
sample path estimate W n (j) of the mean waiting time. The 
mean of these, W n , is an overall estimate of the mean waiting 
time. The medians of the W n (j) are given by the * in the 
box plot. 

Uncorrelated queue (M/M/1) 

N = 1,000 to 10,000 in steps of 1,000 
t = 0.75, 1/E(S) = Service rate = 4.0 
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Box plots derived from 500 replications, j=l,...,500 of the 
sample path estimate W n (j) of the mean waiting time. The 
mean of these, W n , is an overall estimate of the mean waiting 
time. The medians of the W n ( j ) are given by the * in the 
box plot. 

Uncorrelated queue (M/M/1 ) 

N = 11,000 to 20,000 in steps of 1,000 
t = 0.75, 1/E(S) = Service rate = 4.0 
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Figure lc. Box plots derived from 500 replications, j=l,...,500 of the 
sample path estimate Wn(j) of the mean waiting time. The 
mean of these, W n , is an overall estimate of the mean 
waiting time. The medians of the W n (j) are given by the * 
in the box plot. 

Uncorrelated queue (M/M/1 ) 

N = 21,000 to 30,000 in steps of 1,000 
t = 0.75, 1/E(S) = Service rate = 4.0 
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Figure Id. Box plots derived from 500 replications, j=l,...,500 of the 
sample path estimate W n (j) of the mean waiting time. The 
mean of these, W n , is an overall estimate of the mean 
waiting time. The medians of the W n (j) are given by the * 
in the box plot. 

Uncorrelated queue (M/M/1 ) 

N = 31,000 to 40,000 in steps of 1,000 
t = 0.75, 1/E(S) = Service rate = 4.0 
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le. Box plots derived from 500 replications, j=1,...,500 of the 
sample path estimate W n (j) of the mean waiting time. The 
mean of these, W n , is an overall estimate of the mean 
waiting time. The medians of the W n (j) are given by the * 
in the box plot. 

Uncorrelated queue (M/M/1) 

N = 41,000 to 50,000 in steps of 1,000 
t = 0.75, 1/E(S) = Service rate = 4.0 



36 



C. 8500358 
£ 

3 



a 



a 

02 

C 



04 

C4 

X 



- + - 



02 

02 

0 

0 

04 

X 



C2 

0 

02 

0 

c 

0 

X 



0 

02 

0 

02 

03 



0 

0 

0 

0 

02 

0 

C2 

X 



0 

03 

02 

02 

X 



3 


3 




3 






a 




0 






3 


c 




0 




C2 


03 


c 






0 


C2 


0 


03 


0 


C 4 


C2 


X 


C6 


03 


04 




X 


X 


01 



-+- 



-+- -■»- 























1 


1 


j 




1 


1 


• 4 - 






* 




* 








* 




* 




* 




* 




* 




4 




* 










-+ - 
























- 4 - 




- 4 - 



X 



X 



> 


02 


03 


X 




X 


C 4 


C2 


03 


03 


X 


02 


02 


04 


02 


0 


0 


02 


C 


0 


02 


0 


0 




0 


0 


0 


o 


0 










02 


0 


02 






0 


0 


0 






32 


3 






c 


3 


3 











a 



) 

0 

c 

c 

c 

c 



> 

X X 04 

0 02 0 

0 c 

02 C 0 

0 03 02 

0 

0 0 3 

3 

3 3 



C .657*5850 

Figure If. Box plots derived from 500 replications, j=l,...,500 of the 
sample path estimate W n (j) of the mean waiting time. The 
mean of these, W n , is an overall estimate of the mean 
waiting time. The medians of the W n (j) are given by the * 
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Uncorrelated queue (M/M/1) 

N = 51,000 to 60,000 in steps of 1,000 
t = 0.75, 1/E(S) = Service rate = 4.0 
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Figure lg. Box plots derived from 500 replications, j=l,...,500 of the 
sample path estimate W n (j) of the mean waiting time. The 
mean of these, W n , is an overall estimate of the mean 
waiting time. The medians of the W n (j) are given by the * 
in the box plot. 

Uncorrelated queue (M/M/1 ) 

N = 61,000 to 70,000 in steps of 1,000 
t = 0.75, 1/E(S) = Service rate = 4.0 
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Figure 2a. Box plots derived from 500 replications, j=l,...,500 of the 
sample path estimate W n (j) of the mean waiting time. The 
mean of these, W n , is an overall estimate of the mean 
waiting time. The medians of the W n (j) are given by the * 
in the box plots. 

Correlated queue (MDXMD/1) 

N = 1,000 to 10,000 in steps of 1,000 
t = 0.75, 8 = 0.5, p = 0.75. Service rate = 4.0 
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Figure 2b. Box plots derived from 500 replications, j=l,...,500 of the 
sample path estimate W n (j) of the mean waiting time. The 
mean of these, W n , is an overall estimate of the mean 
waiting time. The medians of the W n (j) are given by the * 
in the box plots. 

Correlated queue (MDxMD/1 ) 

N = 11,000 to 20,000 in steps of 1,000 
t = 0.75, 3 = 0.5, p = 0.75. Service rate = 4.0 
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Figure 2c. Box plots derived from 500 replications, j=l,...,500 of the 
sample path estimate W n (j) of the mean waiting time. The 
mean of these, W n , is an overall estimate of the mean 
waiting time. The medians of the W n (j) are given by the * 
in the box plots. 

Correlated queue (MDxMD/1) 

N = 21,000 to 30,000 in steps of 1,000 
t = 0.75, 8 = 0.5, p = 0.75. Service rate = 4.0 
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Figure 2d. Box plots derived from 500 replications, j=l,...,500 of the 
sample path estimate W n (j) of the mean waiting time. The 
mean of these, W n , is an overall estimate of the mean 
waiting time. The medians of the W n (j) are given by the * 
in the box plots. 

Correlated queue (MD^MD/l) 

N = 31,000 to 40,000 in steps of 1,000 
t = 0.75, 8 = 0.5, p = 0.75. Service rate = 4.0 
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Figure 2e. Box plots derived from 500 replications, j=l,...,500 of the 
sample path estimate W n (j) of the mean waiting time. The 
mean of these, W n , is an overall estimate of the mean 
waiting time. The medians of the W n (j) are given by the * 
in the box plots. 

Correlated queue (MDxMD/1) 

N = 41,000 to 50,000 in steps of 1,000 
t = 0.75, S = 0.5, p = 0.75. Service rate = 4.0 
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Figure 2f. Box plots derived from 500 replications, j=l,...,500 of the 
sample path estimate W n (j) of the mean waiting time. The 
mean of these, W n , is an overall estimate of the mean 
waiting time. The medians of the Wn ( j ) are given by the * 
in the box plots. 

Correlated queue (MDxMD/1 ) 

N = 51,000 to 60,000 in steps of 1,000 
t = 0.75, 6 = 0.5, p = 0.75. Service rate = 4.0 
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Figure 2g. Box plots derived from 500 replications, j=l,...,500 of the 
sample path estimate W n (j) of the mean waiting time. The 
mean of these, W n , is an overall estimate of the mean 
waiting time. The medians of the W n (j) are given by the * 
in the box plots. 

Correlated queue (MD^MD/1 ) 

N = 61,000 to 70,000 in steps of 1,000 
t = 0.75, g = 0.5, p = 0.75. Service rate = 4.0 



45 



C. ESTIMATED MEAN WAITING TIMES 



To compare the change of the estimate steady state mean 
waiting time from the M/M/1 queue as p is changed. Table 2a 
shows the value of mean waiting time, W n and the estimate 
standard deviation of those means . The mean service rate 
is fixed at 4.0. 

The mean waiting times of the correlated queue for any 
fixed value of p and t, are plotted in Figures 3a and 3b. 
respectively. These show that the mean waiting time increases 
dramatically as t -> 1, just as for the M/M/1 queue (8 = 1.0), 

y 

for which the mean waiting time increases as l/(l-t). A 
similar drastic increase in the mean waiting time because of 
the increase in cross-correlation, measured by p is apparent 
in Figure 3b, when t is large. Note that when p is close 
to one there is very long term correlation in the queuing 
process; however, p = 1 is not allowed because the system is 
not ergodic. 

The ratio of the estimated mean waiting time, for 
the correlated queue to that for the M/M/1 queue which are 
tabulated in Table 2b, are also plotted in Figure 4a and 
Figure 4b. In particular. Figure 4a shows that the correla- 
tion causes a decrease in mean waiting time for p small, 
and an increase when p is large. The latter effect probably 
occurs because the service times become highly autocorrelated 
as p increases. 

The shape of the plots in Figure 4a suggested the possi- 
bility of fitting the mean waiting time to a function of 
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p/(l-P) and t/(l-t) so as to be able to predict E(W) for 
large values of t. The form of the independent variables 
were . suggested by the fact that a log transform showed 
lnE(W) to go up approximately as ln(l-p) for large p, and 
the fact that for small p, E(W) is smaller than E(W) in 
the M/M/1 case. Also the form t/(l-t) was suggested by the 
M/M/1 theory. 

Only one value of $ was used in the simulations, $ = 0.5, 
so that the results did not reflect this variable. The rough 

fit obtained for t close to one was, for the MDxMD/1 queue 
with parameters t,p,E(S). 

( v . i) e*(w) = iTt [°-5 + ^jrf £ ] E (s) • 

One use for this formula would be to predict E* (W) for 
high t so that a simulation could be started with W q exponen- 
tial with mean E (W ) (the probability of a zero waiting time 

n 

can be neglected for high t) . This reduces the transient 
effect, as we will see later. 

Formula (V.I) predicts that for t = 0.99, p = 0.25 
and E (S) = 0.25, the mean waiting time would be 

E* (W) = 14.43 , t = 0.99, p = 0.25, E(S) = 0.25 

The mean for the M/M/1 queue is given as 
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E (W) 



24.75 



A simulation was performed for this traffic intensity 
(t = .99) , starting with W q = 0; the simulation was found to 
converge to steady-state for n = 280,000 and then (see 
Figures 12a and 12b) 



W 280,000 



14.194 



with 



s.d. ( w 280,000^ 0.1731. 

This is in reasonable agreement. 

The* effect of using this value for E(W q ) in the simulation 
is discussed in Section V.D. . 
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Correlation 



Traffic Intensity, t 





in 

CM 

• 


o 

in 

• 


.75 


.90 


o 

o 

• 


0.073763 

(0.000064) 


0.192998 

(0.000156) 


0.484147 

(0.000501) 


1.258926 

(0.003004) 


. 25 


0.077745 

(0.000730) 


0.210338 

(0.000165) 


0.543177 

(0.000597) 


1.444641 

(0.003373) 


o 

in 

• 


0.083026 

(0.000091) 


0.237792 

(0.000230) 


0.649296 

(0.000899) 


1.798944 

(0.004716) 


in 

• 


0.091539 

(0.000130) 


0.296498 

(0.000440) 


0.926665 

(0.001709) 


2.819429 

(0.009518) 


. 90 


0.101673 

(0.000188) 


0.413239 

(0.000950) 


1.654127 

(0.004754) 


5. 717662 
(0.032587) 


M/M/1 

Queue 


0.083365 

(0.000076) 


0.250211 

(0.000231) 


0.749701 

(0.000983) 


2.254441 

(0.007270) 



Table 2a. Estimated mean W and standard 

deviation of the n estimated mean 

of the average, steady state 

waiting time in the correlated 

queue for various p and t (3 = 0.5 

and_p = 0.25). The estimated s.d. 

of W b is given in brackets, 
n 
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# 



(x Mean Service Time, E ( S ) ) 




Figure 3a. Smoothed, esimated mean waiting time in the correlated 
MDxMD/1 queue as a function of traffic intensity, t. 

In units of mean service time for several values of the 
correlation paramter p. Fixed £ = 0.5. 
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(x Mean Service Time, E(S) ) 




Correlation p 



Figure 3b. Smoothed, estimated mean waiting time in the correlated 

MDxMD/1 queue as_ a^ function of the correlation parameter . 
In units of mean service time for several values of the p 
traffic intensity t. Fixed 6 = 0.5. 
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Correlation 



Traffic Intensity, t 
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.50 


.75 


.90 


. 00 


0.8852 


0.7720 


0.6455 


0.5595 


.25 


0.9329 


0.8414 


0.7242 


0.6421 


o 

If) 
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0.9963 


0.9512 


0.8657 


0.7995 


.75 


1.0985 


1.1860 


1.2356 


1.2531 


. 90 


1.2201 


1.6530 


2.2055 


2.5412 



Table 2b. Ratio of the Estimated mean W n of 
the average waiting time in the 
correlated queue to the known 
average waiting time of the 
uncorrelated (M/M/1) queue (3 = 0.5) 
for various values of p and t. 
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(x E(W) of M/M/1 queue) 




0.00 0.25 0.50 0.75 0.90 1.0 



Traffic Intensity t 

Figure 4a. Smoothed, estimated mean waiting time in the correlated 
MDxMD/1 queue divided by the known mean waiting time in 
the M/M/1 queue as_ a^ function of traf f i c intensity t_. 
Thus the result is in units of the mean waiting time in 
M/M/1 queue. Here the correlation parameter p is the 
parameter. Fixed & = 0.5 
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1 , 

1 . 

0 , 

0 . 



T 



T 




correlation p 



Smoothed, estimated mean waiting time in the correlated 
MDXMD/1 queue divided by the known mean waiting time in 
the M/M/1 queue as a function of the correlation parameter p. 
The result is in units of the mean waiting time in M/M/1 
queue. Here the traffic intensity t is the parameter. 

Fixed 8 = 0.5. 



54 



D. ESTIMATED DISTRIBUTION OF WAITING TIME W 



The simulation program was set up to enable us to 

look at statistics on W and W across the replications 

n n 

at various n's. In addition the file system enabled us 
to look at joing properties of the samples at two different 
n's, and in particular at correlations. The reason for 
this is that a sample of size 500 is rather small when one 
wants to look at the steady state distribution of W n , 
namely W. One could increase the number of replications, 
but this is constrained by the size of the computer. 

However, since the simulation has to be carried out well 
beyond the set in of the steady state for safety sake, 
one might as well pool samples of waiting times along 
the sample path if they are sufficiently far apart to be 
uncorrelated, or approximately so. Then instead of looking 
at a sample size of 500, we may pool waiting times at 10 
steps of n, n^, n^ , . ..,n^g, to make a sample size of 5000. 
Then the waiting time distribution for a particular run can 
be shown as in fig. 5a, 5b; these are actually histograms 
showing many of the estimated sample statistics which 
might be of interest. 

■ In each step of the simulation runs which were performed 
the correlation between waiting times 1000 arrivals apart 
were obtained and found to be negligible, (see Table 3a, 3b) , 
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n = no. of 
arrivals 


6000 


7000 


8000 


9000 


10000 


6000 


1.0 


-0.0118 


0.0135 


0.0002 


-0.0384 


7000 


-0.0118 


1.0 


-0.0450 


0.0403 


-0.1079 


8000 


0.0135 


-0.0450 


1.0 


0.0012 


-0.0088 


9000 


0.0002 


0.0403 


0.0012 


1.0 


-0.0491 


10000 


-0.0384 


-0.1079 


-0.0088 


-0.0491 


1.0 



Table 3a. Correlation Matrix of Waiting Times for 

Different Indices n, in the Correlated Queue, 
t = 0.75, 6 = 0.50, p = 0.75. 
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n = no. of 
arrivals 


6000 


7000 


8000 


9000 


10000 


6000 


1.0 


-0.0371 


0.0280 


0.0105 


-0.0131 


7000 


-0.0371 


1.0 


-0.0413 


0.0779 


-0.1707 


8000 


0.0280 


-0.0413 


1.0 


0.0189 


-0.0241 


9000 


0.0105 


0.0779 


0.0189 


1.0 


-0.0977 


10000 


-0.0131 


-0.1707 


-0.0241 


-0.0977 


1.0 



. Correlation Matrix of Positive Waiting Times 
(No Zeros) for Different Indices n in the 
Correlated Queue. t = .75, 8 = .50, p = 0.75. 



Table 3b 



thus justifying the pooled samples. Of course a more effi- 
cient use of the data would have been to calculate a correlo- 
gram of estimated serial correlations along each (stationary) 
sample path and average over the five hundred replications 
to obtain a more precise estimate of each serial correlation 
together with a variance estimate of that estimate. 

Another reason for looking at the samples W n (j), j = 1, . .., 

500 at different values of n is to assess the convergence to 
a steady state in terms of the whole distribution of W , 
rather than just its mean value. There are a number of ways 
of doing this display and a particularly simple and graphic 
way is to give box-plots for the successive samples on the 
same scale. Unfortunately, this proved to be difficult with 
presently available software;; thus in Figs. 6a-g we give 
box plots of W (j) for n = 1,000(1,000)70,000. The scales 
in these figures are unfortunately not commensurate, but the 
convergence to steady state can be seen. Note that unlike 
the W n (j)'s, the variance of the w n (j)' s are not decreasing 
with n; in fact that should converge to var (W) . If some of 
the distributional detail in the plots is overwhelming, one 
can look solely at the *'s in v the boxes, which are the median 
values in the successive samples. Since these samples are, 
as we have seen, approximately uncorrelated, the *'s could 
have been smoothed. 

The plots in Figures 6a-g are tied to zero by the occurrence 
of zero waiting times (on the average (l-t)% of the sample). 

Thus plots of the positive waiting times are given in Figures 7a-g. 
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Figure 5a. Histogram and sample statistics for combined waiting time 
data in the simulated, uncorrelated M/M/1 queue, i.e. all 
W ( j ) ’ s for j=l , . . . ,500 and n = 51,000(1000)60,000. 
t = 0.75, 1/E(S) = 4.0 
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Figure 5b. Histogram and sample 
data (with 2,293 zero 



time 



statistics for combined waiting 
waiting times removed) in the 
simulated, uncorrelated M/M/1 queue, i.e. all W n (j)'s for 
j=l 500 and n = 51 ,000(1000)60,000. t = 0.75, 1/E(S) = 4.0 
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Figure 6a. Box plots of the sample waiting times W n (j), j=l,...,500 
for values of n = 1,000 to 10,000 in steps of 1,000. 
Correlated queue ; t = 0.75, 6 = 0.5, p = 0.75. 

Figures at top and bottom at left are maximum and minimum 
sample values. 
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Figure 6c. 



Box plots of the sample waiting times W n (j), j=l,...,500 
for values of n = 21,000 to 30,000 in steps of 1,000. 
Correlated queue ; t = 0.75, B = 0.5, p = 0.75. 

Figures at top and bottom at left are maximum and minimum 
sample values. 
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Figure 6d. Box plots of sample waiting times W n (j), j=l,...,500 
for values of n = 31,000 to 40,000 in steps of 1,000. 
Correlated queue ; t = 0.75, 3 = 0.5, p = 0.75. 

Figures at top and bottom at left are maximum and minimum 
sample values. 
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Figure 6e. Box plots of sample waiting times Wp(j), j=l,...,500 
for values of n = 41,000 to 50,000 in steps of 1,000. 
Correlated queue; t = 0.75, 6 = 0.5, p = 0.75. 

Figures at top and bottom at left are maximum and minimum 

sample values. 
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Figure 6f. 



Box plots of sample waiting times W 0 (j), j=l,...,500 
for values of n = 51,000 to 60,000 in steps of 1,000. 
Correlated queue ; t = 0.75, 3 = 0.5, p = 0.75. 

Figures at top and bottom at left are maximum and 
minimum sample values. 
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Figure 6g. Box plots of sample waiting times W (j), j=l,...,500 
for values of n = 61,000 to 70,000 in steps of 1,000. 
Correlated queue ; t = 0.75, B = 0.5, p = 0.75. 

Figures at top and bottom at left are maximum and 
minimum sample values. 
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Figure 7a. Box plots of the sample waiting times W n (j), j=l,. 
wi th zero wai ti nq ti mes removed for values of 
n = 1,000 to 10,000 in steps of 1,000. 

Correlated queue ; t = 0.75, 6 = 0.5, p = 0.75. 
Figures at top and bottom at left are maximum and 
minimum sample values. 
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Figure 7b. Box plots of the sample waiting times W n (j), j = 1 , . 
with zero waiting times removed for values of 
n = 11,000 to 20,000 in steps of 1,000. 

Correlated queue ; t = 0.75, B = 0.5, p = 0.75. 
Figures at top and bottom at left are maximum and 
minimum sample values. 
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Box plots of the sample waiting times W n (j), j=l,. 
with zero waiti ng times removed for values of 
n = 21,000 to 30,000 in steps of 1,000. 

Correlated queue ; t = 0.75, 3 = 0.5, p = 0.75. 
Figures at top and bottom at left are maximum and 
minimum sample values. 
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7d. Box plots of the sample waiting times W (j), j=l,...,500 
with zero waiting ti mes removed for values of 
n = 31,000 to 40,000 in steps of 1,000. 

Correlated queue ; t = 0.75, 3 = 0.5, p = 0.75. 

Figures at top and bottom at left are maximum and 
minimum sample values. 
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7e. Box plots of the sample waiting times W n (j), j=l,...,500 
with zero waiting ti mes removed for val ues of 
n = 41,000 to 50,000 in steps of 1,000. 

Correlated queue ; t - 0.75, 3 = 0.5, p = 0.75. 

Figures at top and bottom at left are maximum and 
minimum sample values. 
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Figure 7f. Box plots of the sample waiting times W n (j), j=l , 
with zero wai ti ng ti mes removed for values of 
n = 51,000 to 60,000 in steps of 1,000. 

Correlated queue ; t = 0.75, 8 = 0.5, p = 0.75. 
Figures at top and bottom at left are maximum 
and minimum sample values. 
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. Box plots of the sample waiting times W n (j), j=l,...,500 
with zero waiting times removed for values of 
n = 61,000 to 70,000 in steps of 1,000. 

Correlated queue ; t = 0.75, 8 = 0.5, p = 0.75. 

Figures at top and bottom at left are maximum 
and minimum sample values. 
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Another question of interest is how the parameters p and 
t affect the waiting time in the MD x MD/1 queue , and how 
these distributions differ from the exponential distribution 
of the positive values of W n for the uncorrelated M/M/1 
queue. There are several graphical techniques for doing this, 
and some of the results are as follows. 

a) Box plots are used to compare the distribution of 
waiting times at 6 = 0.5 for several fixed values of traffic 
intensity t as the correlation parameter p is varied. These 
are shown in fig. 8a-d and complement Figure 3b also for 
fixed values of the correlation parameter p, the variation 
of the distribution of W is shown in Figs. 9a-f. Note that 
in these figures 10 approximately uncorrelated samples in the 
steady state have been combined. 

b) The box plots are not useful for examining departures 

/ 

from exponentiality , and for this purpose the histograms and 
sample statistics from HIST F were utilized. For the grid 
of values of p and t discussed in (V.B) the coefficients of 
variation, skewness and kurtosis were extracted from the 
HISTF outputs (no zeros) , and these are compared to the 
exponential values C(k) = 1.0, y^ = 2.00, y ^ = 6.00 in table 4. 

It is clear that for small p & t the distribution is less 
skewed than the exponential, but as p and t increase the 
distribution becomes positively skewed and is certainly non- 
exponential. Thus for p = 0.9 and t = 0.9 we have C(x) = 1.1325, 
y^ = 2.3289, y^ = 8.6449. One explanation for this is that for 
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Table 4. Deviations of the Estimated Sample Coefficients 
of Variation, Skewness and Kurtosis for the 
Correlated Queue Waiting Times (Without Zeros) 
from the Known (Exponential) Values for the 
M/M/1 Queue. For various p and t values. 
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p small, the service time is decreased by the correlation 
when many people arrive. The distribution is therefore 
tighter than for the M/M/1 case. However as p increases the 
service times become highly correlated through the cross- 
coupling, and this effect, which increases the mean waiting 
time and the skewness of the distribution, is dominant. 

c) Another question of interest is whether for large t 
the heavy traffic approximation holds (Jacobs, 1978) and 
the distribution of W is approximately exponential. This is 
difficult to examine because the time for the queue to reach 
steady state is inordinately high for large t. Thus one case, 
t= 0.995, p = 0.25 was studied and the HISTF outputs are 
shown in Figures 10a & 10b. The plots and the sample statis- 
tics still show a large amount of underdispersion relative 

to an exponential distribution. Thus convergence to an 
exponential distribution as t -»■ 1, if it occurs, must be very 
slow. 

Plots for an extreme case, p = 0, 6=0, are shown 
for various values of t from 0.25 to 0.99 in Figures lla-e. 

d) The question of judging when the transient in the 
simulation has died out is very complex and difficult , just like any 
question of detecting a signal in noise. In particular it is 

to be stressed that the more graphic output one has, the better 
off one is. This is another reason for looking at the W n (j) 
samples as well as the w n (j) samples. 

A 

In Figure 12a is shown for the MD x MD/1 queue 
with t = 0.99, p = 0.25. The judgement was made, and published 
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in Table 1, that equilibrium was reached by the time of the 
40,000th customer and this appears reasonable from the 
figure, which goes out to n = 280,000. The plot of 
is given in Fig. 12b; superposing it on Fig. 12a is a help to 

A 

visualization. It of course has a larger transient than 

because it drags in all the previous annual times , so that the 

/\ 

plot of gives a better picture of the real steady state. 

Of course, as in any real simulation, E(W) is unknown and 
the asymptote in Fig. 12b is a help in showing the convergence. 
Note too that the normality in decreases with n, but that 

/N 

of W does not. 
n 

Note too that since we have determined that W ' s 

n 

1,000 arrivals apart are almost uncorrelated, local smoothing 
could be applied to the W n plot to obtain a better picture. 

Output from a similar simulation for an M/M/1 queue, 
again with t = 0.99 and 500 replications, is shown in Figs. 12c 
St d. The mean of W is known to be E (W) = 24.75. However the 

A 

plot of W n is well above this from approximately n = 60,000 

/\ 

to n = 260,000, and since the correlation between W n values 
is not very extensive, this probably indicates that the 
transient is such that E(W n ) does not go up monotonically to 
E(W), as one might expect intuitively or from the simulation 
out to n = 240,000, but actually rises above E(W) and then 
declines. It may even have a damped oscillation. The point 
is that all the data should be examined- formal tests of the 
similarity of the W n samples at say W^gg ggg and W^g ggg would 
probably have indicated convergence (i.e., similar distributions). 
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Another overlaying the curve of on that of W n 
is a help usually, as is the horizontal asymptotic in Fig. 12d. 
Note that the initial transient in the M/M/1 is longer than 
it is in the MD x MD/1 queue. 

Another question which arises is whether one can 
start the simulation in such a fashion as to reduce the length 
of the transient. An attempt to do this is shown in Figs. 13a-b. 
Here an extreme case, t = 0.995, was taken and the initial 
value W q was not taken to be zero, but an exponential random 
variable with mean 29.0. The value 29.0 comes from the pre- 
diction formula IV. 1. Roughly speaking there appears to be 
faster convergence. The simulation starting from W q = 0 

takes a long time to converge and was too extensive to 
be handled on the computer. 



79 






10.0302544 






32 



a 



a 



7 

9 

9 

C9 



- 4 - 

> 

C.0CC0557 



3 

33 

34 

35 
39 
39 
09 
X 

I 

-+— 

X 



3 

32 

3 

39 

39 

35 

09 

X 

I 

-+ - 

I I 

— * — 

X 



3 

32 
3 2 
34 

34 
37 
36 

35 
35 
05 
X 



I I 
— # — 

X 



.0 



.25 .50 .75 

(CORRELATION ) 



.5 C 



Figure 8a. Correlated queue . Box plots of a sample of waiting times 
Wn-j (j) » j=l , . . . ,500; i=l,...,10; with zero waiting times 
removed, for t = 0.25 and five values of p. For each (t,p) 
pair n-| is chosen to be large enough for the queue to be in 
steady state, and 10 approximately uncorrelated samples at 
n-j ,n 2 ». . . ,n-|Q are combined. 

Here n ] = 31 ,000, 31 ,000, 31 ,000, 31 ,000, 51 ,000 

as p =.0.0, 0.25, 0.50, 0.75, 0.90, respectively. 6 = 0.5. 
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Figure 8b. Correlated queue . Box plots of a sample of waiting times 
Wnj ( J ) » j=l , . . . ,500; i=l,...,10; with zero waiting times 
removed, for t = 0.50 and five values of p. For each (t,p) 
pair n-| is chosen to be large enough for the queue to be in 
steady state, and 10 approximately uncorrelated samples at 
ni ,0^, . . • ,n-|Q are combined. 

Here n ] = 31,000, 41,000, 41,000, 41,000, 61,000 

as p = 0.0, 0.25, 0.50, 0.75, 0.90, respectively . B = 0.5. 
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Figure 8c. Correlated queue . Box plots of a sample of waiting times 
W n ,(j ) > J=1 » • . . ,500; i=l,...,10; with zero waiting times 
removed, for t = 0.75 and five values of p. For each (t,p) 
pair n] is chosen to be large enough for the queue to be in 
steady state, and 10 approximately uncorrelated samples at 
n-| jng,. . . ,n-|Q are combined. 

Here n ] = 41,000, 41,000, 41,000, 61,000, 91,000 

as p = 0.0, 0.25, 0.50, 0.75, 0.90, respectively. 8 = 0.5. 
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Figure 8d. Correlated queue . Box plots of a sample of waiting times 
w nj (377 3=1 >• • • ,500; i=l,...,10; with zero waiting times 
removed, for t = 0.90 and five values of p. For each (t,p) 
pair n-j is chosen to be large enough for the queue to be in 
steady state, and 10 approximately uncorrelated samples at 
ni ,n 2 , . . . ,n-|Q are combined. 

Here n ] = 51,000, 61,000, 61,000, 91,000, 91,000 

as p = 0.0, 0.25, 0.50, 0.75, 0.90, respectively. B = 0.5. 
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Figure 9a. Correlated queue . Box plots of a sample of waiting times 
W n • (j) , j=l ,. . . ,500; i = l,...,10; with zero waiting times 
removed, for p = 0.0 and four values of t. For each (t,p) 
pair n-| is chosen to be large enough for the queue to be 
in steady state, and 10 approximately uncorrelated samples 
at n. jng,. . . ,n.Q are combined, g = 0.5. 



84 



1C. <553716 

a 

3 



a 



a 





a 


5 


3 7 




37 


53 


55 


55 


59 


56 


09 


59 


X 


CS 


1 


X 

I 


1 

- 4 - 


- 4 - 


1 * 1 




- + - 


X 


X 


Q.OCQ0040 




.25 


.50 



52 
5 

5 2 
57 

53 
59. 
59 
59 
55 
C9 
C9 
) 



- + - 




- + - 



X 



.75 

(TRAFFIC INTENSITY) 



5 2 
5 

5 2 
5 2 
5 2 



54 
5 4 
5 3 

53 
5 £ 

54 
5 3 
57 

55 
5 8 
5 5 
55 
5 5 
55 
59 
C 5 
C 5 
09 
C 5 
X 



* 



- 4 - 

I 

X 



.90 



Figure 9b. Correlated queue . Box plots of a sample of waiting times 
• w n-j ( J ) > j=l , . . . ,500; i=1,...,10; with zero waiting times 
removed, for p = 0.25 and four values of t. For each (t,p) 
pair n-| is chosen to be large enough for the queue to be 
in steady state, and 10 approximately uncorrelated samples 
at n-j ,n 2 > . . . ,n-|Q are combined. B = 0.5. 
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Figure 9c. Correlated queue . Box plots of a sample of waiting times 
Wn-j (J ) > J = 1 » • • • .500; i=l,...,10; with zero waiting times 
removed, for p = 0.50 and four values of t. For each (t,p) 
pair n-j is chosen to be large enough for the queue to be 
in steady state, and 10 approximately uncorrelated samples 
at n-j . . . ,n^Q are combined. B = 0.5. 
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Figure 9d. Correlated queue . Box plots of a sample of waiting times 
W n • (j ) , j=l , . . . ,500; i*l,...,10; with zero waiting times 
removed, for p = 0.75 and four values of t. For each (t,p) 
pair n-] is chosen to be large enough for the queue to be 
in steady state, and 10 approximately uncorrelated samples 
at n^ . . ,n^Q are combined. B = 0.5. 
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Figure 9e. Correlated queue . Box plots of a sample of waiting times 
W n . 0 ) » j=l » . • . ,500; i=l,...,10; with zero waiting times 
removed, for p = 0.90 and four values of t. For each (t,p) 
pair n] is chosen to be large enough for the queue to be 
in steady state, and 10 approximately uncorrelated samples 

6 = 0.5. 
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Figure 9f. Correlated queue , p = 0.0, 8 = 0.0 (service time equals 
t times previous inter-arrival time). Box plots of a 
sample of waiting times W n .(j), j=l,...,500; i=l,...,10; 
with zero waiting times removed. Five values of t. For 
each t, n-| is chosen to be large enough for the queue to 
be in steady state, and 10 approximately uncorrelated 
samples at n-j . . . ,n-|g are combined. 
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Figure 10a. 



Histogram and sample statistics for combined waiting-time 
data in the simulated, correlated MDxMD/1 queue; 



t = 0.995, S = 0.5, p = 0.25, 1/E(S) 
i.e. all W n ( j ) ’ s for j=l,...,500 and 



4.0, 

= 91,000(1000)100,000. 
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Figure 10b. Histogram and sample statistics for combined waiting-time 
data <Wi th 18 zero waiting times removed) in the simulated, 
correlated MDxMD/1 queue; t = 0.995. 

1/E(S) = 4.0, i.e. all W (j)’s for ; 



n = 91,000(1000)100,000. 



8 = 0.5, p = 0.25, 
1 , . . . ,500 and 
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11a. Histogram and sample statistics for combined waiting-time 
data (with zero waiting times removed) in the simulated, 
correlated MDxMD/1 queue; t = 0.25, 6 = 0.0, p = 0.0, 
1/E(S) = 4.0, i.e. all W (j)'s for j = 1,...,500 and 
n = 11,000(1000)20,000. n 
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Figure lib. Histogram and sample statistics for combined waiting-time 
data (with zero waiting times removed) in the simulated, 
correlated MDxMD/1 queue; t = 0.50, 3 = 0.0, p = 0.0, 
1/E(S) = 4.0, i.e. all W n (j)'s for j = 1,...,500 and 
n = 21,000(1000)30,000. 
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11c. Histogram and sample statistics for combined waiting-time 
data (with zero waiting times removed) in the simulated, 
correlated MDxMD/1 queue; t = 0.75, 6 = 0.0, p = 0.0, 

1 / E ( S ) = 4.0, i.e. all W n (j)'s for j = 1,...,500 and 
n = 21,000(1000)30,000. " 
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Figure lid. Histogram and sample statistics for combined waiting-time 
data £/ith zero waiting times removed) in the simulated, 
correlated MDxMD/1 queue; t = 0.90, 6 = 0.0, p = 0.0, 
1/E(S) = 4.0, i.e. all W n (j)’s for j = 1,...,500 and 
n = 21,000(1000)30,000. 
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n = 21,000(1000)30,000. 
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CORRELATED QUEUE T=.99 R=.25 B=.50 
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CORRELATED QUEUE T=.99 R=.25 B=.50 
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UNCORRELRTED QUEUE T=.99 
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NUMBER OF ARRIVALS (xl 0^) 

Figure 12c. Estimate W n across 500 replications of the mean waiting time in the M/M/1 queue for 

t = 0.99, and mean service time E(S) = 0.25. The mean appears to fluctuate around the 
known mean (24.75) after 240,000 arrivals. 
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CORRELATED QUEUE T=.995 R=.?5 B=.50 
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VI. VARIANCE REDUCTION IN THE SIMULATION 



A. THE MULTIPLE CONTROL VARIATES METHOD 

Suppose we want to estimate the mean, y of a random 

X 

variable X and suppose V is a zero-mean random variable, 
or can be made so by subtracting from V its known mean. 
Then 



(VI. 1) 



Y = X - aV 



is an unbiased estimator for y and 

x 



(VI. 2) 



a 2 = a 2 + a 2 a 2 - 2a cov[X,V] , 

y x v 



2 2 

where a is the variance of X, a is the variance of 
x v 



V, 



2 2 2 

In order for a to be a minimum given a ,a and cov(X,V], 
y x v 

the optimal choice for a is readily found to be 



(VI. 3) 



a 



* _ Cov [X, V] 



v 



and then, substituting back in (V.III), the minimum variance 
is 



_2 . J2 Cov [X, V] 

y x a 2 

v 



2 2 , 
°x (1 - r x,v> 
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where r is the coefficient of linear correlation between 

X / v 

X and V. Thus the variance of Y, as an estimator of y , 

X 

is always less than that of X if |cov(X,V) | > 1. In a 
simulation X would be the usual estimator of y , and Y would 

X 

be some random variable generated in the simulation which is 
correlated with x. 

2 

Generally, COV (X,V) anda^ are unknown, and we can use in 
its stead an estimate derived from the simulation. This 
idea of reducing the variance of an estimator with a control 
variable is readily extended to a vector V of k control 
variates V^, . . . , with covariance matrix Q. Now let 

Y = X - a ' V 



where a' = (a^, ..., a^,) is a row vector of coefficients, 

and standard results from multivariate analysis show that 



2 2 

a = o +a'Qa-2a'R. 
y x 



Here, Rj = cov[X,Vj]. 



The optimal choice for a is 



a 



* = 



Q 



(assuming, of course that Q is non-singular, for otherwise some 
control variates would be redundant) . Then the minimal value 
of is given by 
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- R'Q 1 R . 



Let M be the joint covariance matrix of the control variates , V 
and X, i . e . , 



For further details see Beja (1968) and Kleijnen (1974).. 

B. SELECTING THE CONTROL VARIABLES 

In order to simulate the mean waiting time in the MDxMD/1 
queue, several control variables derived from the known pro- 
perties of the M/M/1 queue and the service and interarrival 
processes for the MDxMD/1 queue were examined and combined 
to be the multiple control variables. Note that in the simula- 
tion the M/M/1 queue and the MDxMD/1 queue are run with common 
exponential variates, so that the service and interarrival 
times in the M/M/1 queue are the variables from which the 
EARMA type service and interarrival times in the MDxMD/1 
queue are generated. It is well known that in a single-server 
FIFO queue the difference between the cumulative interarrival 




M 



V R Q / 



Then it can be seen that with the optimal 



a 



2 _ det M 

y ~ det Q 
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times and cumulative service times is a good control for W . 

n 

The additional control variables were meant to give any 

greater variance reduction. In the simulation run with 

p = 0.25, t = 0.50 and $ = 0.50, subroutine CONVAR was used 

to compute the correlation between waiting time W , average 

waiting time W n in the MD x MD/1 queue and several quantities 

with known means from uncorrelated queue, i.e., the M/M/1 

waiting time, the M/M/1 average waiting time, the number of 

arrivals since the last regeneration point (W = 0) , sum 

of service time, etc. The controls for W and W are quite 

n n 

different and are examined separately. 

1. Controls for the Waiting Time W n 

The control variables from the uncorrelated queue 

which have been selected as being the most correlated to 

W n are the number of arrivals since the last regeneration 

(V^) , the waiting time (V 2 ) , the average of the previous 100 

interarrival times (V^) , the average of the previous 100 

interarrival times (V^) • The correlation matrix for these 

control variables with W is shown in table 5. 

n 

Since the mean of these four variables are not zero, 
the known expected value must be subtracted from each and 
the estimator is 

Y = X - a' (V - E (V) ) , 



where 
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V = 


(V 1 ,V 2 ,V 3 ,V 4 ) 


a* = 


Pi 
1 — 1 

1 

a 



where 



R = 


cov [X,V] 



and 



X = 


waiting time W n (correlated queue) . 



The values of E (V) used, are as follow: 



E(V 1 ) 


* 

t 

2 

(i-tr 


E(V 2 ) 


t 

= l-t * y s 


e ( v 3 ) 




e(v 4 ) 


= y s • 



Note that even though V^, the sum of the waiting times 
in the MD x MD/1 queue is an average of correlated random 

i k 

See Appendix I. 
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variables, the variables are exponential and their average 
is still p . 

A controlled estimate for was obtained using V^, 
V 2' V 3' V 4 and the resu l ts are shown in Table 6 and Figure 14. 

The 500 replications makes it possible to estimate, as 
in Table 5, the values of the components in a*. In Figure 
14 the second, fourth, ... box plots are for the sample of 
500 controlled W (j" 's, j = 1, ..., 500. There is clearly 
less variability than in the box plots (first, third, etc.) 
for the uncontrolled values W n (j), j = 1/ • .., 500. Successive 
pairs of box plots are for different values of n. 

Table 6 gives, for n = 61,000 (1,000) 70,000 statistics 

of the controlled and uncontrolled W n (j), j = If 500 

samples. The values S(MEAN) are to be compared. Thus for 

n = 70 , 000, S (MEAN) , the estimated standard deviation of W , 

is 0.06375, which is to be compared to the value 0.04632 for 

the controlled sample average. Thus the ratio of the standard 

deviations is 0.04632/0.06375 = 0.726, and the variance 

2 

reduction is (0.726) = 0.528. This is not as much as might 

have been expected from the multiple controls, but it reflects 
the difficulty of getting further control variables which are 
not themselves highly correlated with the previous control 
variables . 

2 . Controls for the Sample-path Average W n 

Different control variables are needed for W n and 
since the. former includes all waiting times up to n and 
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the latter is local. Thus W will only be correlated with 

n 

local controls, while W n will be highly correlated with 

controls which contain the whole history of the sample path. 

Thus 5 control variables from the M/M/1 queue were selected 

to control W time, as follows: 
n 

V number of zero waiting times up until n in the 
M/M/1 queue; 

\? 2 : cumulative interarrival time in the M/M/1 queue; 

V^: cumulative service time in the correlated queue; 

V 4 ? cumulative servie time in the M/M/1 queue; 

V^: mean waiting time, W from the M/M/1 queue. 

The estimated matrix from a simulation run with 500 
replications at p = 0.25, 8 = 0.50, t = 0.50 are shown in 
Table 7 . 

The expected values of V which are used to center V 
are as follows: 

E(V^) = n(l-t) , where n = no. of arrivals; 

n 

E (V- ) = l E(X.) , where E(X.) = mean interarrival 

z i=l 1 1 time of arrivals in 

the correlated queue 

n 

E(V-) = l E (SC ■ ) = ny~, where E(SC.) = mean service time 
J i=l 1 b 1 of arrivals in the 

MD x MD/1 queue; 
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n 



S 



/ 



E(V 4 ) 



n 



l E(S.) 
i=l 1 



where E(S^) = mean service 

time of arrivals 
in the M/M/1 queue 



B(V S ) 



I=t “s 



where = mean service time . 



The results of the simulation run with 500 replications 

and p = 0.75, t = 0.75, 3 = 0.50 are shown in Table 8 and 

Fig. 15. From Table 6 we see that, for n = 70,000, S(MEAN) 

for W n with and without control is, respectively 0.001128 

and 0.001709 and we recall that S (MEAN) is the estimate of 

the standard deviation of W . These values are of course much 

n 

smaller than the corresponding values for W n . Further the 

ratio of the estimated standard deviations is 

0.001128/0.001709 = 0.660035 and the variance reduction is 

2 

(0.660) = 0.4356, a rather healthy value. The variance 

reduction shows up dramatically in Figure 15, which is the 
analog of Figure 14. It is also clear that the controlled. 

W n ' s are symmetric, possibly normally distributed, and that 
they have reached steady state. 
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Control 

Variable 


W 

n 


V l_ 


V 2 


V 3 


V 4 


W 


1.0 


0.1679 


0.4518 


-0.1195 


0.1450 


n 












V 1 


0.1679 


1.0 


0.6657 


-0.2611 


-0.3270 


V 2 


0.4518 


0.6657 


1.0 


-0.2452 


0.0330 


V 3 


-0.1195 


-0.2611 


-0.2452 


1.0 


0.3937 


V 4 


0.1450 


-0. 327 


0.0330 


0.3937 


1.0 



Control Variables 



V, : Number of arrivals since last regeneration in 

M/M/1 queue; 

: Waiting time of M/M/1 queue; 

: Average of 100 previous interarrival times; 

; Average of 100 previous service times in MDxMD/1 
queue . 



Table 5. Correlation Matrix of Waiting Time W n in the 
MDxMD/1 Queue and the Control Variables; 
p = .25, t = .50 and 8 = 0.50. 
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CCFFELATEC CUEUE T = C.75 , e = 0,5 
AVERAGE WAITING TINE 
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Table 8. Estimated Sample Characteristics of Estimated Waiting Times, W (j), j=l, 500 
With and Without Controls at 10 Different Values of n. n 
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Figure 14. Box plots of samples, with and without controls, for estimating the mean waiting time 

in a correlated queue with t = 0.75, 3 = 0.5, p = 0.75. The first, third, ... .nineteenth 
box plots are of W n .(j), j=l,...,500; for m=61,000, n 2 =62,000, ... , ni n =70,000, 
respectively. The second, fourth, ... .twentieth box plots are of the W n .y) with a 
multiple control. Note that the controlled W ni ( j ) may be negative. 1 
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Figure 15. Box plots of samples, with and without controls, for estimating the mean waiting time 
in a correlated queue with t = 0.75, 3 = 0.5, p = 0.75. The first, third, ..., ninet 
box plots are of W n .(j), j=l .... ,500; for ni=61,000, n2=62,000, ... , nin=70,000, 
respectively. The Second, fourth, ..., twentieth box plots are of the w n .(j) with a 
multiple control. 



VII . CONCLUSION 



A. In simulation the queue, the box plots, of the sample 
of waiting times, W n (j), and average waiting time, W (j), 
are helpful in exploring the convergence to steady state. 
Furthermore the box plots of the controlled waiting times, 
compared with the uncontrolled waiting times, are useful 
when examining just how much variance reduction is obtained. 
Note that the variance reduction attained for W n is equivalent 
to cutting down the sample size by at least a half. 

B. An approximation for the mean waiting time for MD x MD/1 

queue in the case of 6 = 0.50 and t + 1 is found to be equal 
to (0.5 + times the mean waiting time of M/M/1 queue. 

Thus the approximate waiting time is 



E (W) = 



1-t 



u s < 0 - 5 + 



C. The multiple control variate technique is useful in re- 

♦ 

ducing the variance of the waiting time and mean waiting time 
when simulating the correlation queue. This is done by simu- 
lating it in parallel with the uncorrelated queue and then 
using the data from the uncorrelated queue to control the 
correlated queue. We found that the effectiveness of the 
variance reduction scheme is better for W n , the sample path 

A 

mean waiting time, than for the waiting time W n * 
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